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Abstract. X-ray spectral imaging provides quantitative imaging of trace elements in 
biological sample with high sensitivity. We propose a novel algorithm to promote the 
signal-to-noise ratio (SNR) of X-ray spectral images that have low photon counts. 
Firstly, we estimate the image data area that belongs to the homogeneous parts 
through confidence interval testing. Then, we apply the Poisson regression through 
its maximum likelihood estimation on this area to estimate the true photon counts 
from the Poisson noise corrupted data. Unlike other denoising methods based on 
regression analysis, we use the bootstrap resampling method to ensure the accuracy of 
regression estimation. Finally, we use a robust local nonparametric regression method 
to estimate the baseline and subsequently subtract it from the X-ray spectral data to 
further improve the SNR of the data. Experiments on several real samples show that 
the proposed method performs better than some state-of-the-art approaches to ensure 
accuracy and precision for quantitative analysis of the different trace elements in a 
standard reference biological sample. 



PACS numbers: 42.30.Va, 87.55.kd, 87.57.cm, 87.64.K-, 87.64.kd 
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1. Introduction 

1.1. Various applications of X-ray spectral imaging 

X-ray spectral imaging has been used for more than half a century to identify and 
quantify the elemental composition of a wide variety of geological, biological, and 
medical target sample (Jenkins et al 1995). Recently, due to the advent of the third 
generation synchrotron radiation facility, X-ray spectral imaging provides quantitative 
imaging of trace elements in biological sample with high sensitivity (sub-mgkg" 1 ) and 
high spatial resolution (sub-um to nm). More and more researchers in the field of 
bio medicine and life science are showing great interest in this technology (Gherase and 
Fleming 2011). In the analysis of diseases such as Parkinson's disease and Alzheimer's 
disease, X-ray spectral images are useful when the quantitative imaging of element 
spatial distribution is needed to study the disease development (Popescu et al 2009, 
Wang et al 2010). Qin et al (2011) used synchrotron radiation X-ray spectral images 
to explore the spatial association of copper in rat aortic media. Furthermore, as an 
interdisciplinary science complementary to genomics and proteomics, a new research 
subject called metallomics has been developed recently and is receiving great attention 
as a new frontier in the investigation of trace elements in biology (Mounicou et al 2009). 
However, the accurate quantitative analysis is badly affected by the Poisson noise and 
baseline errors that are inherent in the X-ray spectral imaging. Especially, denoising 
X-ray spectral images that have low photon counts poses a big challenge in quantitative 
analysis of trace elements in biomedicine, which is also a focus of this study. 

1.2. Typical procedure of X-ray spectral imaging 

Different elements in a sample emit different scattered characteristic X-ray beams of 
many different energies when the sample is scanned and irradiated by the incident 
beams (such as X-ray, electrons) at every scanning location; each of these beams goes 
to the photon counting detector, and the intensities {i.e. photon counts at the detector) 
of these characteristic beams are proportional to the contents of the elements. This 
phenomenon is the very foundation of the analysis based on X-ray spectral imaging. 
Based on this physical law, a typical X-ray energy vs. the intensity spectrum divided 
by thousands of energy channels can be collected. Due to the physical nature of 
characteristic beams, only one or a few elements will be present at a particular energy 
channel of the spectrum when these elements are scanned at particular scanning location. 

Scanning electron microscopy with an energy dispersive X-ray spectrometer (SEM- 
EDS) and energy dispersive micro synchrotron-based X-ray fluorescence (uSXRF) 
imaging are two commonly used methods to study the interactions of trace elements and 
single cells in natural system for the reason that they have relatively high sensitivity 
and high spatial resolution (Twining et al 2003). Under energy- dispersive configuration, 
both methods use the same analytic procedure to quantitatively analyze the spectral 
images. The apparatus of SEM-EDS is cheaper and smaller than that of uSXRF 
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imaging, but the monochromaticity, detection sensitivity and spatial resolution of 
uSXRF imaging are much higher than those of SEM-EDS (Van Grieken and Markowicz 
2002). In this paper, the two data sets used in our experiments are produced by these 
two methods. 

Before using these spectral imaging methods for the quantitative analysis in 
biomedicine, a specific analyte (or standard sample) in the form of either solid or solution 
is placed on the scanning platform to collect the multidimensional X-ray spectral data. 
The acquisition of multidimensional X-ray spectral data is a typical Poisson process 
(Boulanger et al 2010), which makes the raw X-ray spectral data corrupted by Poisson 
noise. Besides, instrument-based systematic errors will also lead to a continuous, slowly 
varying baseline in the acquired spectral data (Twining et al 2003). Fig. 1(a) displays 
such typical raw X-ray spectral data polluted with Poisson noise and systematic baseline. 
Many hardware-based efforts have been made to improve the signal-to-noise ratio (SNR) 
such as the insurance of 90° geometry between the incident and scattered beams (Geraki 
et al 2004). In this paper, we proposed a novel software-based method that can reduce 
Poisson noise and baseline in the X-ray spectral data by means of signal processing. 
The desired effect of our method can be seen from pre-processing procedure in Fig. 1(b) 
with baseline (red line) subtracted from the denoised data (black line). The relatively 
pure spectral data (blue line) in Fig. 1(b) is the output of our proposed method. 

The subsequent procedure is to calculate the intensities (photon counts) of 
characteristic X-ray beams of different elements from the spectrum. In this step, the 
signals of different elements need to be separated and then integrated (Fig. 1(c)). 
The original data (black dash line) has been separated as the sum of Ca (red peaks), 
Fe (green peaks), Cu (blue peaks) and Zn (purple peaks). Methods such as iterative 
least-squares fitting of a mathematical model combined with Monte Carlo simulations 
(Bekemans et al 2003), baseline-corrected spectra fitting to a summed exponentially 
modified Gaussian (EMG) peak model with a sigmoidal baseline (Twining et al 2003) 
are typically used in the separation of different elements' characteristic peaks. After 
the separation procedure, each element's peak areas are integrated to be equal to the 
number of photon counts for each element, which is then used for the quantitative 
mapping of element spatial distribution to disease development (Fig. 1(d)). The whole 
data processing procedure are displayed in Fig. 1. 

Recently, a multi-platform open source software called PyMCA for the analysis 
of energy- dispersive X-ray fluorescence spectra has been developed (Sole et al 2006). 
This software has combined many well performed algorithms in it, which can be used 
to calculate the intensities of characteristic X-ray beams (as displayed in Fig. 1(c)) in 
our standard biological samples. In next step (Fig. 1(d)), we demonstrate the use of 
these characteristic photon intensities so that we can obtain the quantitative amounts 
of different trace elements in standard biological samples (Wang et al 2010). 
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Figure 1. The typical processing procedure for X-ray spectral image data: (a) 
Simulated raw X-ray spectral data which are corrupted with Poisson noise and 
baseline; (b) Using preprocessing methods to reduce Poisson noise and baseline, baseline 
(red line) is subtracted from the denoised data (black line) and the clean signal 
(blue line) is obtained (here the baseline is lowered manually in order to give a clear 
show);(c) Separating different characteristic peaks and using these peaks to calculate 
the characteristic X-ray intensities of different elements, here we simulated the KL peak 
of Ca, KL and KM peaks of Fe, Cu, and Zn; the preprocessed data (black dash line) 
is raised manually in order to give a clear displayed) the calculated X-ray intensities 
can be used in the subsequent analysis such as mapping certain elements in biological 
samples; 



1.3. Review of Poisson denoising and baseline removal 

As has been mentioned above, we will describe a new signal processing method that 
deals with Poisson noise and baseline in the X-ray spectral image data. Our method is 
generally applicable to X-ray spectral images that have low photon counts and therefore 
pose a big challenge in quantitative analysis of trace elements in biomedicine. 

There is an extensive literature on Poisson denoising methods which can be 
generally divided into three classes. The first use multiscale analysis technique (Zhang 
et al 2008, Luisier and Blu 2008, Wang 2007, Spring and Clegg 2009) such as wavelet 
analysis. After the noisy signal being decomposed into noise and the useful signal by 
wavelet transform, the inversely transformed signals will be free from noise. Since the 
Poisson statistics are generally more difficult to be tackled than the Gaussian ones, 
the variance stabilizing transform is integrated into the multiscale analysis framework 
to transform the noise model from Poisson to Gaussian (Anscombe 1948, Spring and 
Clegg 2009, Makitalo and Foi 2011, Zhang et al 2008, Palakkal and Prabhu 2012). In 
general, these algorithms usually require a prior knowledge of noise to set up appropriate 
parameters. The second class of methods estimate the true photon counts directly 
through statistical means, such as Bayesian inference combined with multiscale analysis 
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(Timmerman and Nowak 1999, Kolaczyk 1999, Lefkimmiatis et al 2009), hypothesis 
testing (Kolaczyk E D 2000), maximum likelihood estimation and regression analysis. 
In addition, there are variational approaches for Poisson denoising (Le et al 2007, Chan 
and Chen 2007, Bonettini and Ruggiero 2011, Zhou and Li 2013). Regression analysis 
methods (Boulanger et al 2008, Kervrann and Trubuil 2004) have a good performance 
when they deal with low photon count data. However, the statistical methods often 
give bad results when the sample size is small. Some resampling algorithms (Haynor 
and Woods 1989) have been introduced to overcome this small sample problem. Among 
these algorithms, bootstrap method (Dahlbom 2002) allows estimation of the sampling 
distribution of almost any statistic using only very simple methods. Considering the 
multidimensional X-ray spectral data are usually low photon count data containing a 
very high level of Poisson noise, we propose Poisson regression with bootstrap resampling 
methods to remove the Poisson noise in the multidimensional X-ray spectral data. To 
the best of our knowledge, we are the first to propose the bootstrap Poisson regression 
methods for X-ray spectral image denoising. 

In order to further improve the performance of X-ray spectral imaging, we should 
also remove the baseline in the X-ray spectra. The baseline is caused by systematic errors 
such as the non-linear response of the detector (Van Grieken and Markowicz 2002). The 
continuous and low varying baseline in X-ray spectra can be treated as superposition 
on the original data (Ruckstuhl et al 2001). Various methods in literature have been 
proposed to estimate the baseline of a spectrum. Liland et al (2010) give an overview 
of the baseline correction for multivariate calibration of spectra. Here we choose the 
robust nonparametric regression methods based on the algorithm proposed by Ruckstuhl 
et al (2001) for baseline removal since this method is both effective and robust. Being 
different from the original method, our algorithm uses a new regression kernel for the 
flexibility of second-order robust regression. 

The remainder of this paper is organized as follows. Section 2 will explain the 
method in details. Section 2.2 introduces the bootstrap Poisson regression method. 
Section 2.3 outlines the robust regression methods for baseline removal. Section 3 
gives the experimental results of our method in comparison with other state-of-the-art 
methods using standard samples including biological samples. In Section 4 we briefly 
summarize our method and future research directions. 

2. Methods 

2.1. Data acquisition and data structure 

In this work, we use alloy wire data and two standard biological samples (bovine liver 
NIST 1577a and pig liver GBW 08551) as real samples in our multidimensional X- 
ray spectral imaging experiments. Both wire sample and biological samples are placed 
on an acquisition platform with micron spot of incident beam focused on the each 
scanning location in the samples. An energy dispersive X-ray spectrometer counts the 
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characteristic photons emitted from each scanning sample point. Then the incident 
beam spot will move to the next scanning location for continuous data acquisition. The 
scanning time of each scanning point is the same. Based on the above description, 
the underlying structure of the acquired multidimensional X-ray spectral data can be 
taken as a 2D-1D structure. The 2D part of the multidimensional data refers to the two 
dimensional images that are formed from the total scanning points of the same energy 
channel while the ID part refers to the whole spectrum for all energy channels at a 
single scanning point. 

2.2. Poisson denoising 

Photon counting X-ray spectral image data are typical Poisson distributed data. That 
is to say the photon count data Yj (i = 1, ■ ■ • , N, with N being the total number of 
scanning points in each image of the 2D part of the data) follow a distribution with 

density function A») = - ± yr\ — , where Aj is the desired noise free photon counts, 
which can be estimated as mathematical expectation of Yj. However, it is impossible 
to calculate the expectation of f(Yf, Aj) (A*) through Yi itself. One method to solve the 
problem is to find more data from the same distribution (/(Y^; Aj)). As for X-ray spectral 
data, these identically distributed data can be found through the following analysis. 

In modern high-resolution X-ray spectral imaging, an incident beam is scanned as a 
nanoscale spot in a raster pattern across the sample's surface to make the scanning points 
get fairly close to each other to offer sufficient details of the sample. The nanoscale spot is 
so small that the neighbouring scanning points around an interest point can be assumed 
to belong to a homogeneous region of the sample. In the biomedical imaging, the 
homogeneous regions are not equally extended to all directions from a point of interest, 
so that these irregular homogeneous regions usually introduce local discontinuities at 
some directions in the sample. Therefore, we need to check which neighbouring points 
belong to the same homogeneous part as the current scanning point of interest. A 
proper size of neighborhood should be determined so that there is sufficient number 
of homogeneous points that are chosen to recover the noise free photon counts Aj. 
Furthermore, the bootstrap resampling method is used to produce sufficient candidates 
for good estimation of Aj. With these spatial-domain local homogeneity assumption 
and sampling strategies, we have local photon count data that are independent and 
identically distributed (IID) for the following Poisson regression analysis. 

By analyzing the above-mentioned spatial-domain features of the X-ray spectral 
data, we can use Poisson regression analysis to estimate the Aj from these IID data. 
Poisson regression assumes the response variable lj has a Poisson distribution, and 
assumes the logarithm of its expected value (Aj) can be modeled by a linear combination 
of unknown parameters. We suppose that the total number of homogeneous data in the 
local Poisson regression area is m, and these homogeneous data that belong to the same 
distribution as that of Yi in the neighbourhood of lj are marked as Yy (Yj is included 
in Yij,j = !,•••, m). Therefore, the logarithmic forms of the expectation Aj for the Yjj 
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can be fitted with a linear model function: 

log(Xi) = aiXj + k + e 4 ; A; « Aj = e a ^ +b ',j = 1, • • • , m (1) 

where is an auxiliary explanatory variable for each point of Yjj, A, is the estimator 
of the true expectation Aj, a, and &j are unknown parameters to be calculated, £j is 
the white noise with a fixed variance and zero mean. The explanatory variable Xj has 
no physics meaning but acts as a mathematic auxiliary tool for Poisson regression. To 
determine the explanatory variable Xj, the Yjj are sorted increasingly or decreasingly so 
that the explanatory variable Xj can simply be the sorting index. Therefore, the noise 
free photon counts Aj are approximated as: 

Xi = E{Yi) « e aiX i' +bi (2) 

where xy is the corresponding explanatory variable to Yj. 

With the above-mentioned general scheme in mind, we should first choose the 
proper irregular homogeneous regions so that the Yj's neighbouring photon count data 
have the same distribution function as that of Yj. Here we use the significance test 
to solve this problem. The acquired photon counts Yj are assumed to have a Poisson 
distribution. The data Y i in the Yj's neighbourhood should share the same cumulative 
distribution function: 

m = e~* £ TT (3) 

k=0 

With equation (3) we can calculate the lower and upper bound photon counts of the 
confidence interval of level 95%. The data Y- that have the photon counts outside the 
confidence interval will be taken as the data samples having the different distributions 
from that of Yj. 

After choosing the homogeneous data that are candidate for choosing Yjj in Poisson 
regression analysis of Yj, the performance of estimator Aj in equation (1) is then 
dependent on the Yij's size (or bandwidth) m and can vary at each point of the photon 
count data sequence according to image contents. Small bandwidth will make local 
regression analysis sensitive to noises and outliers, while large bandwidth will create 
a large approximation error in local regression. In order to optimally estimate the 
bandwidth m, we analyze the performance of the estimator and consider the usual local 
L 2 risk (Boulanger et al 2008) defined as: 

K(\ i ,\ i ) = E[(\ i -\ i ) 2 ] (4) 

where Aj is the unknown expectation. The local risk lZ(\i, Aj) is defined at each point 
and then differs from usual global performance measures that integrate errors on the 
whole images. This local risk of the candidates selected from the significance test reaches 
its minimum at each point. Choosing a new larger m that does not increase the number 
of candidates means that the local discontinuity appears; the previous smaller m is then 
considered as the optimal size. Otherwise, the local risk should be minimized to choose 
the optimal m as follows. 
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Boulanger et al (2008) have given in detail the solution of the minimization 
problem described in equation (4), which designs a sequence of increasing bandwidth: 
M = {m^(xj),n E [0, N); m^ n ~ l \xj) < m^ n '(xj)}. And then this sequence is used to 
detect the optimal bandwidth m* for the local smoothing: 

m*(xj) = sup {n' < n : \\ n ^ — Xf 1 ^ < dv n i(x)} (5) 

m(")(i 3 )eA/ 

where $ = 2-*/2 is a positive constant, v 2 ,(x) is the variance of the data Yi in the 
regression context. Boulanger et al (2008) have proved that the m* in the sequence of 
M that satisfies equation (5) will be the optimal bandwidth m that minimizes the risk 
described in equation (4). Typically, choosing eight or more homogeneous data points 
in the local regression area will ensure satisfactory estimation of Aj. 

In order to calculate the values of aj and 6, in the equation (1), we use the principle 
of maximum likelihood estimation to compute the set of parameters (a*, bi) that make 
the following log-likelihood function value as large as possible: 

m 

/(a,, k) = y. J">-':i + h) - e aiXj+bi - log{Y id \) (6) 
i=i 

Unfortunately, directly computing Oj and bi is difficult since equation (6) has no closed- 
form solution. An iterative weighted least squares method (Davison and Hinkley 1997) 
can be used to estimate a\ and bi. At each iteration an adjusted responses vector 
Zi = (• • • , Zij, • • •) is regressed on the Xj with elements z it j being expressed as 



z %i = w i + ( y i,j -e Wi ) * — (7) 

Wj 

where the weight Wj is given by the estimators a« and bi (corresponding to and bi) 

Wj = hiXj + b\ (8) 
The results of each iteration given in the form of matrix is 

= (X T WX)~ 1 X T X^ (9) 



0>i 

bi 



where X is the matrix of Xj, W is the diagonal matrix of weights Wj (i.e. W[s, t] = Wj, 
when s — t; W[s, t] = 0, when s ^ t). 

So far, the estimators aj and bi are ready to be substituted into the equation (2) 
for estimating the noise free photon counts Aj. However, the accuracy of statistical 
estimation depends on the sample size. In practice, one needs as many samples as 
possible to ensure high estimation accuracy in statistical analysis. However, it is often 
not easy to obtain many samples. Efron (1979) has introduced bootstrap resampling 
methods to deal with this problem. Here we also use the bootstrap resampling methods 
to increase the accuracy of statistic analysis. 

To apply bootstrap method to enhance the accuracy of estimating Aj, we need to 
generate new sample data Y*j from the original data Specifically, we generate Y*^ 
from the estimated expectation Aj by using following equation: 

Y^ = Xi + >i /2 *e}J = l,---,m (10) 
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where £*,••• £^ are adjustment parameters that determine the performance of the whole 
resampling method. It has been demonstrated that the residuals of Poisson regression 
can be used as these adjustment parameters for bootstrap resampling methods (Davison 
and Hinkley 1997). In order to deduce the residuals of Poisson regression, we first 
introduce the H matrix that is derived from equations (7)-(9): 

H = X(X T WX)~ 1 X T W (11) 

With this definition of matrix H, the standardized Pearson residuals can be written as 

Y — A 

r Pj = — — hvS'-? = 1 >"""> m ( 12 ) 

{A,(i-/gr 

where hj is the jth diagonal element of the hat matrix H. 

The standardized Pearson residuals are further expressed as mean-adjusted Pearson 
residuals by rpj — fp, where fp is the mean of rp.y These mean-adjusted Pearson 
residuals have all the qualities that bootstrap resampling method needs. With the 
adjustment parameters £*,••• e* m being sampled from these mean-adjusted, standardized 
Pearson residuals according to the bootstrap resampling rule, the new sample data Y*j 
can be obtained as additional data to implement Poisson regression. Thus the desired 
expected value Aj (noise free photon counts in Equation (2)) can be estimated with 
higher accuracy. Theoretically the number of new samples generated from bootstrap 
methods needs to be infinity. In practice, 300 or more samples can ensure good results 
with pleasant accuracy. 



2.3. Baseline drift removal 

Baseline drift in X-ray spectral images is caused by hardware-based systematic errors 
such as the non-linear response of the detector. Based on the feature of slowly-varying 
local continuous baseline, it is convenient to remove baseline drift from the spectrum of 
each scanning point by using a robust local regression method (Ruckstuhl et al 2001). 

After the Poisson regression analysis, the spectral data of one scanning point can 
be modeled as follows: 

V(c k ) = g(c k ) + s(c fe ) +e h ,k = l,---K (13) 

where V(c k ) is the processed data by Poisson regression, g(c k ) is the baseline, s(c k ) is 
the desired signal at the energy channel c k , and e k represents the measurement errors 
with zero mean and variance £, K is the total number of points in the whole spectrum 
of a single scanning point. 

In order to separate the three components in equation (13), a locally weighted 
scatter plot smoothing method (Cleveland 1979) can be used. With the data V(c k ) at a 
energy channel c k , we suppose the bandwidth of the local regression context around c k 
is n. In the n defined local regression context, the Poisson regression processed data t k+i 
within this context can be defined with the following modified equation (from equation 
(13)): 

t k+i = g{c k+i ) + E k+i \ E k+i = s(c k+i ) + e k+i \ i = -n/2, ■ ■ ■ , n/2 (14) 
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(15) 



where g(c k+ i) is the baseline function that has enough smoothness and E k+i is the sum 
of signal s(c k+i ) and error e k+i . The baseline function g(c) can be approximated by 
using second-order Taylor's formula at point c k : 

g(c) = g(c k ) + g'(c k )(c - c k ) + g"{c k )(c - c k ) 2 + 0((c - c k ^ 2 

g(c) « g(c) = Pq + Pi{c - c fc ) + /3 2 (c - c k 

where g{c) is the baseline estimated by the regression model. The parameter vector 
P(c k ) = [p , fa, P2Y can be calculated by incorporating a weight scheme into the local 
least-squares problem to decreases the influence of data points in proportion to their 
distance from c k . That is, 

n/2 

P(c k ) = argmin J2 K (^T ■){**+< _ \Po + /M c fc+i - c fc) 

^ i=-n./2 11 

+(3 2 (c k+i -c k ) 2 ]} (16) 

where i^[(cfc+j — c k )/h] is a unimodal symmetric nonnegative weight function that is 
zero outside the c k s neighbourhood, which is defined by c k ± h with h being half of 
the bandwidth of the regression context. The estimation performance is not greatly 
dependent on the choice of the weight function K (Ruckstuhl et al 2001). Here we 
choose a logarithmic function to reduce the influence of the quadratic part in equation 
(15), so that the weight of data in the energy channels far from the current channel will 
drop quickly to have little effect on the estimation while the data in the close energy 
channels will have high weight: 

K{u) = 1 - log[(e - l)u + 1] (17) 

Equation (16) is used to estimate f3(c k ) initially, next we use the residuals of this 
estimation to assign robustness weight w r (c k+ i) to each point, such that the points 
with large residuals will receive small robustness weights. The baseline curve g(c) is 
then refined by performing a weighted least-squares fit, according to 

n/2 

P(c k ) = argmin w r {c k+i )K{^- ){t k+i - [p + Pi{c k+i - c k ) 

i=—n/2 

+P 2 (c k+i -c k ) 2 }} (18) 

This fit is repeated iteratively to converge with the weights w r (c k+ i) always being 
calculated from the previous iteration. On the choice of various w r (c k+ i), we use Tukey's 
bisquare weights 

w r (c k+i ) = {max[l - (r k+i /b) 2 , 0]} (19) 

where r k+i = [t k+i — g(c k+i )]/a. The a is estimated using the standardized median of 
absolute values of the residuals 

a = median(\t k+ i — ^(cfc+j)|)/0.6745 (20) 

The whole baseline estimation is summarized as follows: 
(a) Equation (16) is used to compute g(c k ) initially; 
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(b) Use equation (19) to calculate the robustness weights w r (ck+i); 

(c) Use equation (18) to compute a new fitted value g(c.k). 

The iterative steps (b) and (c) are repeated until \g n +i(ck) — g n {ck)\ < 0.01, where 
g n (ck) and <? n +i(cfc) are the results of step (c) at the nth and the (n + l)th iteration, 
respectively. Usually ten iterations will be enough. 

The local regression bandwidth n is the only parameter that needs to be set. For 
small size n, the robust local regression estimator is more likely to estimate g(ck) + s(ck) 
than g{ck). To avoid such a failure, a value n being 2.5 times the full width of the widest 
peak in the spectra is recommended. 

3. Results and Analysis 

In this section, we firstly use alloy wire experiment data to test the performance of our 
method by comparing it with the four state-of-the-art algorithms. Then we perform 
quantitative analysis on two standard biological samples (bovine liver NIST 1577a and 
pig liver GBW 08551) by using our algorithm and the two best algorithms from the four 
state-of-the-art algorithms. 

3.1. Alloy wire experiment 

3.1.1. Experimental data Alloy wire sample is described in figure 2 which consists of 
a series of six types of wires embedded in an epoxy matrix, with the wire alloys being 
composed from a pallet of six different elements. This alloy sample is imaged with SEM- 
EDS. Typical data acquisition conditions for this type of spectral image are described 
in detail by Kotula et al (2003) and more details about this sample can be found in 
Keenan (2007). Figure 2(a) shows a standard SEM-EDS image of this sample together 
with the composition and component abundances. The image dimensions are 128 x 128 
pixels, and a complete 1024-channel spectrum is collected at each pixel. Figure 2(b) 
shows a typical single-pixel spectrum for the Cu/Mn/Ni wire. The discrete nature of 
the data is clearly evident, and the SNR is sufficiently low so that the presence of Ni 
can not be clearly discriminated from the background. 

3.1.2. Results and analysis In order to test the bootstrap Poisson regression with 
robust nonparametric regression baseline removal (BPR-RR) method, we introduce 
other four state-of-the-art methods for comparison purpose. The first method uses 
traditional Anscombe transform combined with Wiener filter (ATW) to remove the 
Poisson noise and local medians (LM) algorithm to remove the baseline artifacts 
(Friedrichs 1995). All parameters of this ATW-LM method are set as follows: the 
local filtering area of Wiener filter is 3 x 3; the window width of local medians is 200. 
The second method is designed for denoising mixed Poisson and Gaussian noise (MPG) 
(Zhang et al 2008) and removing baseline by asymmetric least squares (ALS) methods 
(Eilers 2004). The parameters of MPG- ALS method are selected as: Poisson noise's 
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Figure 2. (a) The different positions of the wires with the different wire compositions 
and component concentrations: (a) 100% Ni (b) 36% Ni, 64% Fe (c) 70% Cu, 30% Zn 
(d) 16% Cr, 84% Fe (e) 13% Mn, 4% Ni, 83% Cu (f) 100% Cu; ((a)-(f) arc referred 
to the rows of dots in figure 2(a)). (b) A single-pixel spectrum from the Cu/Mn/Ni 
Wire. 



weight alpha is set to 1; the mean and standard deviation of Gaussian component are 
and 5; the number of iterations is set to 10. The SURE-LET (Luisier and Blu 2008) and 
BLS-GSM (Wang 2007) method are also introduced to denoise multidimensional X-ray 
spectral data. We set the parameters of these two methods as recommended in the 
original papers. The parameters of our proposed BPR-RR method are: The number 
of bootstrap samples is set to 300; the bandwidth of the local regression context for 
baseline estimation is set to 250. With those parameters all the methods mentioned 
above achieve their best performances. 

Since the whole spectra data contain 1024 spectral images, it is not possible to 
display all of them in this paper. Here we choose the typical spectral images at KL2 
line energy channel to show the denoising and baseline drift removal performance with 
the colour scale representing different numbers of photon counts. We will also give 
the corresponding ID plots showing single-pixel spectra for different elements. From 
the visual inspection aspect, the performance can be evaluated in terms of removing 
the noises and baseline drift at background regions and preserving the original photon 
counts at target regions. 
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Figure 3. The X-ray spectral images at the energy channel of Mn's KL2 Line Energy 
(5.89 KeV) for (a) original data, (b) BLS-GSM method, (c) ATW-LM method, (d) 
MPG-ALS method, (e) SURE-LET method, (f) BPR-RR method. 
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Figure 4. Typical single-pixel spectra of Mn/Ni/Cu alloy wire (line e in figure 2(a), 
the line displayed in figure 3) for (a) original data, (b) BLS-GSM method, (c) ATW- 
LM method, (d) MPG-ALS method, (e) SURE-LET method, (f) BPR-RR method 
(vertical axis: photon counts; horizontal axis: cV). 
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Figure 5. The X-ray spectral images at the energy channel of Cr's KL2 Line Energy 
(5.41 KeV) for (a) original data, (b) BLS-GSM method, (c) ATW-LM method, (d) 
MPG-ALS method, (e) SURE-LET method, (f) BPR-RR method. 
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Figure 6. Typical single-pixel spectra of Cr/Fe alloy wire (line d in figure 2(a), the 
line displayed in figure 5) for (a) original data, (b) BLS-GSM method, (c) ATW-LM 
method, (d) MPG-ALS method, (c) SURE-LET method, (f ) BPR-RR method (vertical 
axis: photon counts; horizontal axis: cV). 
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Figure 7. The X-ray spectral images at the energy channel of Zn's KL2 Line Energy 
(8.62 KcV) for (a) original data, (b) BLS-GSM method, (c) ATW-LM method, (d) 
MPG-ALS method, (e) SURE-LET method, (f) BPR-RR method. 



(a) (b) (c) 




Figure 8. Typical single-pixel spectra of Cu/Zn alloy wire (line c in figure 2(a), the 
line displayed in figure 7) for (a) original data, (b) BLS-GSM method, (c) ATW-LM 
method, (d) MPG-ALS method, (e) SURE-LET method, (f ) BPR-RR method (vertical 
axis: photon counts; horizontal axis: eV). 
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Figure 9. The X-ray spectral images at the energy channel of Fe's KL2 Line Energy 
(6.39 KeV) for (a) original data, (b) BLS-GSM method, (c) ATW-LM method, (d) 
MPG-ALS method, (e) SURE-LET method, (f) BPR-RR method. 

The low-concentration elements Mn, Cr and Zn have very weak photon counts and 
low SNRs in their original data (see figure 3(a), 5(a) and 7(a)), which are more apparent 
in the corresponding single-pixel spectra (figure 4(a), 6(a) and 8(a)). After denoising, the 
MPG-ALS method have changed the original photon counts of these low-concentration 
elements at target regions (see figure 3(d), 5(d) and 7(d)), and the spectral waveforms 
have distorted (figure 4(d), 6(d), and 8(d)). The BLS-GSM method smoothes the raw 
spectral data too much as the photon counts at target regions have been homogenized 
(see figure 3(b), 5(b) and 7(b)) while preserves relatively strong noises in the spectra 
(figure 4(b), 6(b) and 8(b)). The denoised data by SURE-LET method resemble the 
original photon count data very much for the positions of elements Mn, Cr and Zn, 
while the denoising performance is not desirable at background and other component 
elements' positions (see figure 3(e), 5(e) and 7(e)). SURE-LET method also has the 
worst baseline drift removing performance compared with the other methods (figure 
4(e), 6(e) and 8(e)). 

Among all the five methods, BPR-RR and ATW-LM methods achieve the best 
performance in terms of both denoising and preserving the original photon counts. 
Figure 3(c)(f), 5(c)(f) and 7(c)(f) show the satisfying denoising effect of BPR-RR and 
ATW-LM methods for element Mn, Cr and Zn simultaneously with the high fidelity to 
original photon counts. However, we can see that there are less noise left on the spectra 
of BPR-RR method than the spectra of ATW-LM method (figure 4(c) (f), 6(c) (f) and 
8(c) (f)). Although the BPR-RR method shows a decrease in photon count (intensity) 
compared to the raw data (figure 7(a) and (f)), this performance is acceptable since the 
raw data are contaminated with Poisson noise, which means the raw maximum does not 
always represent the real maximum photon counts. In single-pixel X-ray spectra, the 
KM peak of Cu (the peak at 890 on figure 8. (a)(f)) and the KL peak of Zn (the peak at 
862 on figure 8. (a)(f)) are so close that the peak of Cu will cause an effect of baseline 
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Figure 10. The signal-to-noise ratios estimated from the X-ray spectral images of 
each element's KL2 line energy. The elements from the left to right are Cr (5.41 KeV), 
Mn (5.89 KeV), Fe (6.39 KeV), Ni (7.46 KeV), Cu (8.03 KeV) and Zn (8.62 KeV). 



drift on the peak of Zn, which means the energy channel 862 contains two elements (Cu 
and Zn) rather than one element (Zn). For this reason, BPR-RR method has partly 
removed the baseline drift effect of Cu and thus causes the loss of denoised photon 
counts of Zn (figure 7(f)). As for element Fe with high SNR, all the above methods 
have similar good results except SURE-LET method (see figure 9). The corresponding 
ID plots of single-pixel spectra of the lower line in figure 9 is actually displayed in figure 
6. 

Although the noiseless data is unavailable, we decided to estimate the SNR over 
the whole data. Keenan (2007) has provided a way to calculate the SNR of Poisson 
data. With this method, the SNR of the original X-ray spectral data matrix D can be 
calculated as 

gjujj^ sum °f eigenvalues of D T D 
sum of all data elements in D 
However, the equation (21) can not be used to calculate the SNR of spectra data 
preprocessed by above-mentioned methods because the preprocessed data have lost their 
Poisson nature. All of the preprocessed data can be seen as true data corrupted with 
Gaussian white noise. Since no wire signal is present in the background areas (black areas 
in figure 2(a)), data of these clear background areas can be taken as noise to calculate 
the variance of Gaussian noise in the preprocessed data. However, this variance of the 
clear background area cannot be used as the variance of the ensemble image unless 
the noise process is ergodic. Actually, the acquisition of whole multidimensional X- 
ray photon count data and inherent Poisson noise follows a typical Poisson process 
(Boulanger et al 2010) with ergodic properties (Wolff 1981). Furthermore, our space- 
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independent denoising algorithm will not change the ergodicity. Therefore, we use the 
variance of the background the noise variance of the ensemble image. By means 

of standard image processing methods such as threshold, it is easy to pick up the data 
in the background area. The signal mean is easily calculated by using the data on the 
areas of wire dots. Finally, the SNR of preprocessed data is computed as a ratio of the 
signal mean to the noise standard deviation (the square root of noise variance). 

Figure 10 shows the SNRs of the images in the KL2 line energy channels of the 
six component elements. Figure 10 display that component elements with small SNRs 
still have the relative small SNRs in the processed data. The MPG-ALS method and 
BPR-RR method can enhance the SNRs of low-concentration elements (Cr, Mn and 
Zn) to almost the same level as the SNRs of high-concentration elements (Ni, Fe and 
Cu). The ATW-LM method, SURE-LET method and BPR-RR method have similar 
high computed SNRs while our method has the highest computed SNRs among these 
methods. (All the elements' line energy data are referred to the database of the National 
Institute of Standards and Technology of U. S.) 

3.2. Standard biological sample experiment 

3.2.1. Experimental data The two commercially available standard biological samples, 
bovine liver (NIST 1577a) and pig liver (GBW 08551) with known element 
concentrations, are imaged with uSXRF at the beamline BL15U at Shanghai 
Synchrotron Radiation Facility (Shanghai, China). These two samples were separately 
pressed into a ~ 5mgcm -2 tablet, and then sandwiched between double Mylar films. 
The image dimensions of bovine liver and pig liver are 11 x 11 and 10 x 10, respectively. 
Both samples are scanned with a complete 2048-channel spectrum at each pixel. More 
details about the beamline station and the preparation of samples can be found in Wang 
et al (2010). 

3.2.2. Results and analysis Based on the results of section 3.1, we choose ATW-LM 
method, SURE-LET method and our BPR-RR method to preprocess the two standard 
biological samples. However, SURE-LET method failed to process the spectra data 
of these samples because the quantity of these two samples data are too small to be 
processed by SURE-LET method. All the parameters of the other methods were set as 
in section 3.1. 

We further use the method described in Marco et al (1999) to implement 
quantitative analysis of these two biological samples. In uSXRF imaging, the varied 
current intensity of X-ray and the differences in thickness and density of thin biological 
samples can result in considerable errors and undesirable precision. Thus, for enhancing 
quantification precision, an internal reference in the test specimen must be used to 
demonstrate the relationship between the fluorescent intensity of the analyte and a 
signal from an internal reference. As proposed in Marco et al (1999), the Compton peak 
of uSXRF spectra is possible to be used as an internal standard for trace element (e.g., 
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Ca, Fe, Cu and Zn) quantification in organic matter due to the Compton scattering 
being theoretically related to the mass of the sample. The commercially available 
standard reference materials, having closely matched matrix with the analyte and known 
element concentration, are used to compute the sensitivity R{ C of element % relative to 
the Compton peak using the following equation: 

Ric = y~tt- ( 22 ) 

l cr \^i T 

where Jj r is the fluorescent intensity of each element % in the standard reference sample, 
I cr is the intensity for the Compton scattering peak of the standard reference sample, 
Ci r is the known element concentrations of the standard reference sample. 

With the relative ratio Ri c of the sensitivity of each element i, the concentration 
Ci a of each element in the analyte can be calculated by the intensity of each element % 
and Compton scattering peak area of the analyte, i.e. 

^ = T^T (23) 

J-ca -t<4c 

where ij a is the fluorescent intensity of each element % in the analyte, I ca is the intensity 
for the Compton scattering peak of the analyte. 

Based on the above-described method, we choose the bovine liver as the standard 
reference sample with the matrix-matched pig liver being treated as the standard analyte 
to valuate our proposed method. According to the element compositions in these two 
standard biological samples, we choose the elements Ca, Fe, Cu and Zn to perform the 
quantitative calculations. Firstly, we use three different bovine liver data, i.e., the BPR- 
RR, ATW-LM method preprocessed bovine liver data and the raw bovine liver data, 
to calculate the three different sets of the intensities for each element Ij r and Compton 
peak I cr by peak area integration (figure 1(c)), and then compute the three different 
sets of sensitivities R{ C of element % relative to the Compton peak according to the 
equation (22) (see table 1). Although the same bovine liver sample making specimen 
matrices and element contents constant in the three preprocessed data, the different 
preprocessing methods with (or without) different Poisson denoising and baseline drift 
removal methods can affect the relative sensitivity i?j C . 

Secondly, we use the BPR-RR, ATW-LM method preprocessed pig liver data and 
the raw pig liver data to calculate the three different sets of intensities for each element 
Ii a and the Compton peak I ca . Then by using the Ri C calculated from the matrix- 
matched bovine liver data, the concentrations of Ca, Fe, Cu and Zn in the standard pig 
liver data were calculated according to the equation (23) and in comparison with the 
certified values in table 2. 

All the intensities for the selected elements and Compton scattering peak are 
calculated by PyMCA. Although PyMCA itself has preprocessing methods, here we 
only use this software to separate and integrate different elements' characteristic peaks 
and Compton peaks in the raw liver data, BPR-RR and ATW-LM method preprocessed 
liver data. 
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Table 1. The ratios Ri C obtained from the raw bovine liver (NIST 1577a) data and 
the data preprocessed by the ATW-LM and BPR-RR methods 



Ric x 10 2 




raw data 


ATW-LM 


BPR-RR 


Ca 


0.75 ±0.04 


0.69 ±0.04 


0.77 ±0.05 


Fe 


1.47 ±0.15 


1.50 ±0.15 


1.59 ±0.16 


Cu 


2.73 ±0.12 


2.68 ±0.12 


2.86 ±0.13 


Zn 


3.36 ±0.22 


3.47 ±0.23 


3.62 ±0.24 



Table 2. The concentrations of elements (ugg x ) in pig liver data (GBW 08551) 
calculated by using the different Ri C from the raw and the two preprocessed bovine 
liver data (NIST 1577a) at table 1. 





Concentrations of elements: ugg 


-l 




Certified value 


raw data 


ATW-LM 


BPR-RR 


Ca 


197 ±14 


256 ± 15 


235 ± 14 


195 ±11 


Fe 


1050 ± 40 


1144 ± 118 


1062 ±97 


1055 ± 109 


Cu 


17.2 ± 1.0 


16.2 ±0.7 


16.5 ±0.7 


17.4 ±0.8 


Zn 


172 ±8 


199 ±13 


167 ±8 


178 ± 12 



Table 2 shows that Poisson denoising and baseline drift removal methods do have 
positive effect on the accuracy and precision for quantitative analysis of the different 
trace elements. Among the three different preprocessing methods, BPR-RR method 
has produced best quantitative analysis of element concentrations that are in good 
agreement with the certified values. The BPR-RR method has obtained high accuracy 
and high precision for quantitative analysis of the trace element Ca and Cu. Although 
the quantitative analysis of Zn is not as good as those of Ca and Cu, the BPR-RR 
method still has good accuracy and good precision for quantitative analysis of trace 
element Zn. The possible cause of the decrease of accuracy and precision for quantitative 
analysis of Zn is already mentioned in the section 3.1.2. Both bovine liver and pig liver 
contain the elements of Cu and Zn. The peaks of the spectra of Cu and Zn are so 
close that even an excellent tool such as PyMCA is not able to separate these two 
elements in the spectra thoroughly and precisely. As for the element Fe, the BPR-RR 
method has achieved good accuracy but poor precision for quantitative analysis due 
to the calculated standard deviation uncertainty (±109) being significantly larger than 
the certified standard deviation uncertainty (±40). One explanation for this is that the 
concentration of element Fe in pig liver (GBW 08551) is very high to introduce large 
photon counts in the spectrum. Large photon counts assume an approximately normal 
distribution rather than Poisson distribution (Haight 1967). Therefore, the method that 
is designed to deal with the low photon count data {e.g. Cu, Ca and Zn) will not have 
the excellent denoising effect for the high photon count data of element Fe at pig liver. 



Reducing Poisson noise and baseline drift 



20 



4. Conclusion 

In this paper, we have introduced a bootstrap Poisson regression and robust 
nonparametric regression method to reduce Poisson noise and baseline drift in the 
multidimensional X-ray spectral data. The proposed method is very effective to improve 
the SNRs of the raw X-ray spectral data. The comparison with other competing methods 
shows that the BPR-RR method offers performance better than some state-of-the-art 
approaches. By using two standard biological samples and applying Compton peak 
standardization in uSXRF quantitative imaging, the BRP-RR preprocessing method can 
ensure satisfying accuracy and precision for quantitative analysis of the trace elements 
in biological samples. The concentrations of Ca, Fe, Cu and Zn in the standard reference 
material (GBW 08551, pig liver) determined by BRP-RR preprocessing and subsequent 
quantitative imaging method were in good agreement with the certified values. This 
work can be extended along several directions in the future. First, the bootstrap samples 
can be sampled through nonparametric ways to converge more quickly to the statistic 
features of the original raw data. Second, Poisson regression can be computed within 
the Bayesian framework, where a priori information about the expectation of photon 
count data can be introduced to improve its performance. Third, our algorithm cannot 
achieve satisfactory performance in the case of extreme local discontinuities (or local 
nonhomogeneity) at all directions from the current scanning point of interest in the 
sample, though, such case is a rare occurrence in modern high-resolution biomedical 
X-ray spectral imaging. At last, the Poisson denoising and baseline removal method 
are automatically adaptive not only to the low-concentration elements but also to the 
high-concentration elements that produce high photon count data (such as the element 
Fe in liver data), while our proposed method is apt to deal with the low photon count 
data. 
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